Effect of felzartamab anti-CD38 treatment on the molecular phenotype of antibody-mediated rejection in kidney transplant biopsies.

1 Analysing the effect of treatment on molecular scores using aligned-rank transformed analysis of variance

1.1 INTRODUCTION

This analysis assesses how felzartamab therapy affects molecular sores measured with gene expression data from biopsies from kidney allografts taken from patients with antibody mediated rejection. The assumptions of normality and heterogeneity of variance differ across molecular scores, thus we avail to alinged-rank transformation to carry out a non-parametric ANOVA. This is powerful in this context because it allows us to test for the interactive effect of treatment (i.e., placebo vs felzartamab) and time (i.e., change in scores from biopsy to biopsy). Aligned-rank transformation also allows for us to specify a mixed-model to handle repeated measurements (i.e., biopsies) across patients.

References:

• Wobbrock J, F.L., Gergle D, Higgins J The Aligned Rank Transform for Nonparametric Factorial Analyses Using Only ANOVA Procedures. In Proceedings of the ACM Conference on Human Factors in Computing Systems (CHI ’11) 143–146. doi:10.1145/1978942.1978963, https://depts.washington.edu/acelab/proj/art/.(2011).

• Kay M, E.L., Higgins J, Wobbrock J. ARTool: Aligned Rank Transform for Nonparametric Factorial ANOVAs. doi:10.5281/zenodo.594511, R package version 0.11.1, https://github.com/mjskay/ARTool. (2021).

• Elkin, L.A., Kay, M., Higgins, J.J. & Wobbrock, J.O. An Aligned Rank Transform Procedure for Multifactor Contrast Tests. in The 34th Annual ACM Symposium on User Interface Software and Technology 754–768 (Association for Computing Machinery, Virtual Event, USA, 2021).

1.2 DATA WRANGLING

First we need to load the necessary R packages and data.

# HOUSEKEEPING ####
# Load necessary libraries
library(tidyverse)
library(broom)
library(rstatix)
library(flextable)
library(officer)
library(emmeans)
library(multcomp)
library(ARTool)
library(rcompanion)
library(Biobase)
library(ggpubr) 
library(patchwork) 
library(ggprism) 

# Custom operators
"%nin%" <- function(a, b) match(a, b, nomatch = 0) == 0

# Suppress dplyr reframe info
options(dplyr.reframe.inform = FALSE)

# Load data
load("C:/R/CD38-effect-of-treatment/natmed/data/cd38_3Sept24.RData")
load("C:/R/CD38-effect-of-treatment/natmed/results/flextable_artANOVA.RData")
# load plot data
load("C:/R/CD38-effect-of-treatment/natmed/results/plots_artANOVA.RData")

# Seed for reproducibility
set.seed(42)

Now we define the categories of molecular scores for easier organization downstream.

# Feature categories
vars_cfDNA <- c("cfDNA_cpml")
vars_abmr <- c("ABMRpm", "ggt0", "ptcgt0", "DSAST", "NKB")
vars_tcmr <- c("TCMRt", "tgt1", "igt1", "QCAT", "TCB")
vars_injury <- c("IRRAT30", "IRITD3", "IRITD5", "cigt1", "ctgt1")
vars_parenchyma <- c("KT1", "KT2")
vars_macrophage <- c("AMAT1", "QCMAT")

# DEFINE VARIABLES TO ASSESS ####
vars <- c(vars_cfDNA, vars_abmr, vars_tcmr, vars_macrophage, vars_injury, vars_parenchyma)

Now we process the data so we can iteratively carry out the analyses on each score.

# DEFINE THE DATA ####
data <- set %>%
    pData()  %>% 
    dplyr::select(Center, Patient, Felzartamab, Group, Followup, Felzartamab_Group, Felzartamab_Followup, all_of(vars)) %>%
    dplyr::filter(Patient %nin% c(15, 18)) %>%
    left_join(., summarise(., sample_pairs = n(), .by = Felzartamab_Followup), by = "Felzartamab_Followup") %>%
    relocate(sample_pairs, .after = Felzartamab_Followup) %>%
    arrange(Felzartamab, Patient, Group)


# WRANGLE THE PHENOTYPE DATA ####
data_00 <- data %>%
    expand_grid(category = c("cfDNA", "ABMR", "TCMR", "macrophage", "injury", "parenchyma")) %>%
    nest(.by = category) %>%
    mutate(
        features = map(
            category,
            function(category) {
                if (category == "cfDNA") {
                    vars_cfDNA
                } else if (category == "ABMR") {
                    vars_abmr
                } else if (category == "TCMR") {
                    vars_tcmr
                } else if (category == "macrophage") {
                    vars_macrophage
                } else if (category == "injury") {
                    vars_injury
                } else if (category == "parenchyma") {
                    vars_parenchyma
                } 
            }
        ),
        data = pmap(
            list(features, data),
            function(features, data) {
                data %>%
                    dplyr::select(
                        Center, Patient, Felzartamab, Group, Followup, Felzartamab_Group, Felzartamab_Followup, sample_pairs,
                        all_of(features)
                    )
            }
        )
    )


# WRANGLE THE DATA FOR UNIVARIATE TESTS ####
data_01 <- data_00 %>%
    mutate(
        data_univariate = pmap(
            list(data, features),
            function(data, features) {
                data %>%
                    pivot_longer(
                        cols = all_of(features),
                        names_to = "variable",
                        values_to = "value"
                    ) %>%
                    nest(.by = variable) %>%
                    mutate(
                        variable = variable %>%
                            factor(
                                levels = vars
                            ),
                        annotation = case_when(
                            variable %in% vars_cfDNA ~ "cfDNA",
                            variable %in% vars_tcmr ~ "TCMR-related",
                            variable %in% vars_abmr ~ "ABMR-related",
                            variable %in% vars_macrophage ~ "macrophage-related",
                            variable %in% vars_injury ~ "injury-related",
                            variable %in% vars_parenchyma ~ "parenchyma-related",
                            TRUE ~ " "
                        ) %>%
                            factor(
                                levels = c(
                                    "cfDNA",
                                    "ABMR-related",
                                    "TCMR-related",
                                    "macrophage-related",
                                    "injury-related",
                                    "parenchyma-related",
                                    "archetypes",
                                    "rejection PC",
                                    "injury PC"
                                )
                            ),
                        score = case_when(
                            variable == "cfDNA_cpml" ~ "Donor-derived cell-free DNA (dd-cfDNA, cp/mL)",
                            variable == "TCMRt" ~ "TCMR classifier (TCMRProb)",
                            variable == "TCB" ~ "T cell burden (TCB)",
                            variable == "tgt1" ~ "Tubulitis classifier (t>1Prob)",
                            variable == "igt1" ~ "Interstitial infiltrate classifier (i>1Prob)",
                            variable == "QCAT" ~ "Cytotoxic T cell-associated (QCAT)",
                            variable == "ABMRpm" ~ "ABMR classifier (ABMRProb)",
                            variable == "DSAST" ~ "DSA-selective transcripts (DSAST)",
                            variable == "NKB" ~ "NK cell burden (NKB)",
                            variable == "ggt0" ~ "Glomerulitis classifier (g>0Prob)",
                            variable == "ptcgt0" ~ "Peritubular capillaritis classifier (ptc>0Prob)",
                            variable == "AMAT1" ~ "Alternatively activated macrophage (AMAT1)",
                            variable == "QCMAT" ~ "Constitutive macrophage (QCMAT)",
                            variable == "IRITD3" ~ "Injury-repair induced, day 3 (IRITD3)",
                            variable == "IRITD5" ~ "Injury-repair induced, day 5 (IRITD5)",
                            variable == "IRRAT30" ~ "Injury-repair associated (IRRAT30)",
                            variable == "cigt1" ~ "Fibrosis classifier (ci>1Prob)",
                            variable == "ctgt1" ~ "Atrophy classifier (ct>1Prob)",
                            variable == "KT1" ~ "Kidney parenchymal (KT1)",
                            variable == "KT2" ~ "Kidney parenchymal - no solute carriers (KT2)"
                        ), .before = 1
                    ) %>%
                    arrange(annotation, variable)
            }
        )
    ) %>%
    dplyr::select(category, data_univariate) %>%
    unnest(data_univariate) %>%
    mutate(n_cat = score %>% unique() %>% length(), .by = category, .after = variable)

This is what the data look like after we have organized them. Each row of the dataframe contains the data for each score, the category it belongs to, its category annotation, its long-format name, and the name of the variable in the data.

data_01 %>% print(n="all")

## # A tibble: 20 x 6
##    category   annotation         score         variable n_cat data    
##    <chr>      <fct>              <chr>         <fct>    <int> <list>  
##  1 cfDNA      cfDNA              Donor-derive~ cfDNA_c~     1 <tibble>
##  2 ABMR       ABMR-related       ABMR classif~ ABMRpm       5 <tibble>
##  3 ABMR       ABMR-related       Glomerulitis~ ggt0         5 <tibble>
##  4 ABMR       ABMR-related       Peritubular ~ ptcgt0       5 <tibble>
##  5 ABMR       ABMR-related       DSA-selectiv~ DSAST        5 <tibble>
##  6 ABMR       ABMR-related       NK cell burd~ NKB          5 <tibble>
##  7 TCMR       TCMR-related       TCMR classif~ TCMRt        5 <tibble>
##  8 TCMR       TCMR-related       Tubulitis cl~ tgt1         5 <tibble>
##  9 TCMR       TCMR-related       Interstitial~ igt1         5 <tibble>
## 10 TCMR       TCMR-related       Cytotoxic T ~ QCAT         5 <tibble>
## 11 TCMR       TCMR-related       T cell burde~ TCB          5 <tibble>
## 12 macrophage macrophage-related Alternativel~ AMAT1        2 <tibble>
## 13 macrophage macrophage-related Constitutive~ QCMAT        2 <tibble>
## 14 injury     injury-related     Injury-repai~ IRRAT30      5 <tibble>
## 15 injury     injury-related     Injury-repai~ IRITD3       5 <tibble>
## 16 injury     injury-related     Injury-repai~ IRITD5       5 <tibble>
## 17 injury     injury-related     Fibrosis cla~ cigt1        5 <tibble>
## 18 injury     injury-related     Atrophy clas~ ctgt1        5 <tibble>
## 19 parenchyma parenchyma-related Kidney paren~ KT1          2 <tibble>
## 20 parenchyma parenchyma-related Kidney paren~ KT2          2 <tibble>

1.3 CALCULATE MEDIAN SCORES AND MEDIAN EFFECT OF TREATMENT

At this stage we can summarize the median values for each score by treatment and follow-up.

(Note that the median scores are used for general interpretation of treatment effects because aligned-ranks, the data being assessed for inference regarding the effects of treatment, are not easily interpretable.)

# UNIVARIATE MEDIANS ####
data_02 <- data_01 %>%
    mutate(
        medians = map(
            data,
            function(data) {
                data %>%
                    reframe(
                        median = value %>% median(),
                        IQR = value %>% IQR(),
                        sample_pairs = n(),
                        .by = c(Followup, Felzartamab, Felzartamab_Followup)
                    ) %>%
                    arrange(Felzartamab_Followup)
            }
        ),
        medians_delta = map(
            medians,
            function(medians) {
                medians %>%
                    reframe(
                        median_delta = combn(median, 2, diff) %>% as.numeric(),
                        IQR_delta = combn(IQR, 2, mean) %>% as.numeric(),
                        sample_pairs = sample_pairs,
                        .by = c(Felzartamab)
                    ) %>%
                    mutate(
                        Followup_pairwise = rep(c("Baseline - Week24", "Baseline - Week52", "Week24 - Week52"), 2) %>%
                            factor(levels = c("Baseline - Week24", "Baseline - Week52", "Week24 - Week52")),
                        .before = sample_pairs
                    ) %>%
                    mutate(
                        median_delta_delta = combn(median_delta, 2, diff) %>% as.numeric(),
                        IQR_delta_delta = combn(IQR_delta, 2, mean) %>% as.numeric(),
                        .by = c(Followup_pairwise),
                        .before = sample_pairs
                    )
            }
        )
    )

Let’s take a look at what the median summaries look like for the ABMRprob score.

data_02 %>%
    dplyr::filter(variable == "ABMRpm") %>%
    pull(medians) %>%
    pluck(1) %>%
    flextable::flextable()

Followup	Felzartamab	Felzartamab_Followup	median	IQR	sample_pairs
Baseline	Placebo	Baseline_Placebo	0.6695	0.35525	10
Week24	Placebo	Week24_Placebo	0.7740	0.39025	10
Week52	Placebo	Week52_Placebo	0.5185	0.49825	10
Baseline	Felzartamab	Baseline_Felzartamab	0.7360	0.30775	10
Week24	Felzartamab	Week24_Felzartamab	0.1675	0.35700	10
Week52	Felzartamab	Week52_Felzartamab	0.7025	0.53875	10

Now we can look at the change in median ABMRprob scores in the placebo and felzartamab arms, as well as the differnce in those medians (i.e., ΔΔ median - this is our summarized treatment effect).

data_02 %>%
    dplyr::filter(variable == "ABMRpm") %>%
    pull(medians_delta) %>%
    pluck(1) %>%
    flextable::flextable()

Felzartamab	median_delta	IQR_delta	Followup_pairwise	median_delta_delta	IQR_delta_delta	sample_pairs
Placebo	0.1045	0.372750	Baseline - Week24	-0.6730	0.3525625	10
Placebo	-0.1510	0.426750	Baseline - Week52	0.1175	0.4250000	10
Placebo	-0.2555	0.444250	Week24 - Week52	0.7905	0.4460625	10
Felzartamab	-0.5685	0.332375	Baseline - Week24	-0.6730	0.3525625	10
Felzartamab	-0.0335	0.423250	Baseline - Week52	0.1175	0.4250000	10
Felzartamab	0.5350	0.447875	Week24 - Week52	0.7905	0.4460625	10

1.4 CARRY OUT ALIGNED-RANK TRANSFORM ANOVA

Now we carry out the analysis. This includes testing the specific contrasts for the interaction of treatment and time (i.e., our effect of treatment).

# UNIVARIATE NONPARAMETRIC TESTS ####
artANOVA <- data_02 %>%
    mutate(
        art = map(data, art, formula = value ~ Followup * Felzartamab + (1 | Patient)),
        art_aov = map(art, anova, type = "II"),
        art_aov_tidy = map(art_aov, tidy) %>% suppressWarnings(),
        art_aov_contrast = map(
            art,
            art.con,
            formula = "Followup:Felzartamab", adjust = "fdr", method = "pairwise", interaction = TRUE, response = "art"
        ),
        art_aov_contrast_tidy = map(art_aov_contrast, tidy),
        art_aov_contrast_cld = map(
            art_aov_contrast_tidy,
            function(art_aov_contrast_tidy) {
                art_aov_contrast_tidy %>%
                    as.data.frame() %>%
                    cldList(adj.p.value ~ Followup:Felzartamab, data = .)
            }
        )
    )

1.5 CREATE A SUMMARY TABLE OF ANOVA RESULTS FOR EACH SCORE

# CREATE FLEXTABLE SUMMARY OF ART MODELS ####
title_art <- paste("Table i. Non-parametric ANOVA (ART) of molecular scores in biopsies from treated vs untreated patients")
artANOVA %>%
    dplyr::select(annotation, n_cat, score, art_aov) %>%
    unnest(everything()) %>%
    dplyr::rename(p.value = `Pr(>F)`) %>%
    mutate(FDR = p.value %>% p.adjust(method = "fdr"), .by = annotation) %>%
    dplyr::select(annotation, score, Term, `F`, p.value, FDR) %>%
    flextable::flextable() %>%
    flextable::add_header_row(values = rep(title_art, ncol_keys(.))) %>%
    flextable::merge_h(part = "header") %>%
    flextable::merge_v(j = 1:2) %>%
    flextable::fontsize(size = 8, part = "all") %>%
    flextable::align(align = "center", part = "all") %>%
    flextable::bg(bg = "white", part = "all") %>%
    flextable::bg(i = ~ FDR < 0.05 & Term == "Followup:Felzartamab", j = 3:6, bg = "#fbff00") %>%
    flextable::colformat_double(j = 2:3, digits = 2) %>%
    flextable::colformat_double(j = 4:ncol_keys(.), digits = 3) %>%
    flextable::border_remove() %>%
    flextable::bold(part = "header") %>%
    flextable::padding(padding = 0, part = "all") %>%
    flextable::border(border = officer::fp_border(), part = "all")  %>%
    flextable::autofit()

Table i. Non-parametric ANOVA (ART) of molecular scores in biopsies from treated vs untreated patients
annotation	score	Term	F	p.value	FDR
cfDNA	Donor-derived cell-free DNA (dd-cfDNA, cp/mL)	Followup	4.418	0.019	0.029
		Felzartamab	1.598	0.222	0.222
		Followup:Felzartamab	10.396	0.000	0.001
ABMR-related	ABMR classifier (ABMRProb)	Followup	6.659	0.003	0.007
		Felzartamab	1.620	0.219	0.235
		Followup:Felzartamab	9.605	0.000	0.001
	Glomerulitis classifier (g>0Prob)	Followup	11.768	0.000	0.001
		Felzartamab	1.511	0.235	0.235
		Followup:Felzartamab	10.743	0.000	0.001
	Peritubular capillaritis classifier (ptc>0Prob)	Followup	14.188	0.000	0.000
		Felzartamab	3.178	0.092	0.125
		Followup:Felzartamab	7.749	0.002	0.004
	DSA-selective transcripts (DSAST)	Followup	10.508	0.000	0.001
		Felzartamab	1.747	0.203	0.234
		Followup:Felzartamab	4.179	0.023	0.035
	NK cell burden (NKB)	Followup	6.512	0.004	0.007
		Felzartamab	1.956	0.179	0.224
		Followup:Felzartamab	5.114	0.011	0.018
TCMR-related	TCMR classifier (TCMRProb)	Followup	3.555	0.039	0.165
		Felzartamab	1.417	0.249	0.468
		Followup:Felzartamab	1.124	0.336	0.504
	Tubulitis classifier (t>1Prob)	Followup	1.016	0.372	0.506
		Felzartamab	1.853	0.190	0.408
		Followup:Felzartamab	0.065	0.937	0.968
	Interstitial infiltrate classifier (i>1Prob)	Followup	3.414	0.044	0.165
		Felzartamab	0.997	0.331	0.504
		Followup:Felzartamab	0.032	0.968	0.968
	Cytotoxic T cell-associated (QCAT)	Followup	5.693	0.007	0.107
		Felzartamab	4.118	0.057	0.172
		Followup:Felzartamab	0.927	0.405	0.506
	T cell burden (TCB)	Followup	4.360	0.020	0.151
		Felzartamab	2.868	0.108	0.269
		Followup:Felzartamab	0.301	0.742	0.856
macrophage-related	Alternatively activated macrophage (AMAT1)	Followup	1.549	0.226	0.304
		Felzartamab	1.047	0.320	0.320
		Followup:Felzartamab	2.406	0.105	0.209
	Constitutive macrophage (QCMAT)	Followup	3.378	0.045	0.209
		Felzartamab	3.100	0.095	0.209
		Followup:Felzartamab	1.426	0.253	0.304
injury-related	Injury-repair associated (IRRAT30)	Followup	0.055	0.946	0.946
		Felzartamab	3.842	0.066	0.223
		Followup:Felzartamab	1.915	0.162	0.347
	Injury-repair induced, day 3 (IRITD3)	Followup	0.174	0.841	0.946
		Felzartamab	3.268	0.087	0.223
		Followup:Felzartamab	3.125	0.056	0.223
	Injury-repair induced, day 5 (IRITD5)	Followup	0.073	0.930	0.946
		Felzartamab	1.088	0.311	0.583
		Followup:Felzartamab	2.586	0.089	0.223
	Fibrosis classifier (ci>1Prob)	Followup	0.360	0.700	0.876
		Felzartamab	3.482	0.078	0.223
		Followup:Felzartamab	0.938	0.401	0.668
	Atrophy classifier (ct>1Prob)	Followup	0.568	0.572	0.780
		Felzartamab	3.297	0.086	0.223
		Followup:Felzartamab	0.602	0.553	0.780
parenchyma-related	Kidney parenchymal (KT1)	Followup	0.968	0.389	0.518
		Felzartamab	1.628	0.218	0.518
		Followup:Felzartamab	0.605	0.552	0.552
	Kidney parenchymal - no solute carriers (KT2)	Followup	0.859	0.432	0.518
		Felzartamab	0.911	0.352	0.518
		Followup:Felzartamab	1.038	0.365	0.518

1.6 SUMMARISE THE EFFECT OF TREATMENT FOR EACH SCORE ACROSS ALL TIME PERIODS

Here we present the effect of treatment (ΔΔ) for each score, along with the effect of time (Δ) for placebo and felzartmab arms.

flextable_pairwise

Table 1. Effect of felzartamab on molecular scores
Annotation	Score	Baseline - Week 24				Week 24 - Week 52				Baseline - Week 52
		Δ Placebo (N=10)	Δ Felzartamab (N=10)	ΔΔ	ΔΔ FDR	Δ Placebo (N=10)	Δ Felzartamab (N=10)	ΔΔ	ΔΔ FDR	Δ Placebo (N=10)	Δ Felzartamab (N=10)	ΔΔ	ΔΔ FDR
cfDNA	Donor-derived cell-free DNA (dd-cfDNA, cp/mL)	0.50 (76.12)	-55.00 (30)	-55.50 (53.06)	2e-04	-3.50 (69.12)	14.50 (25.5)	18.00 (47.31)	0.027	-3.00 (51)	-40.50 (45.75)	-37.50 (48.38)	0.045
ABMR-related	ABMR classifier (ABMRProb)	0.10 (0.37)	-0.57 (0.33)	-0.67 (0.35)	1e-03	-0.26 (0.44)	0.54 (0.45)	0.79 (0.45)	1e-03	-0.15 (0.43)	-0.03 (0.42)	0.12 (0.43)	0.709
	Glomerulitis classifier (g>0Prob)	-0.02 (0.23)	-0.31 (0.35)	-0.29 (0.29)	3e-03	-0.19 (0.34)	0.31 (0.41)	0.50 (0.37)	2e-04	-0.22 (0.3)	0.00 (0.38)	0.21 (0.34)	0.277
	Peritubular capillaritis classifier (ptc>0Prob)	-0.12 (0.25)	-0.55 (0.27)	-0.43 (0.26)	3e-03	-0.10 (0.37)	0.37 (0.31)	0.47 (0.34)	3e-03	-0.22 (0.3)	-0.18 (0.35)	0.04 (0.33)	0.977
	DSA-selective transcripts (DSAST)	-0.09 (0.22)	-0.37 (0.24)	-0.28 (0.23)	0.026	-0.16 (0.32)	0.16 (0.33)	0.32 (0.33)	0.026	-0.24 (0.33)	-0.20 (0.25)	0.04 (0.29)	0.988
	NK cell burden (NKB)	0.02 (0.22)	-0.78 (0.49)	-0.80 (0.36)	0.014	-0.19 (0.4)	0.68 (0.59)	0.87 (0.49)	0.014	-0.17 (0.46)	-0.10 (0.62)	0.07 (0.54)	0.972
TCMR-related	TCMR classifier (TCMRProb)	-0.01 (0.02)	0.00 (0)	0.01 (0.01)	0.358	0.00 (0.02)	0.00 (0)	0.00 (0.01)	0.856	-0.01 (0.02)	0.00 (0)	0.01 (0.01)	0.358
	Tubulitis classifier (t>1Prob)	-0.01 (0.06)	-0.02 (0.04)	-0.01 (0.05)	1.000	0.00 (0.04)	0.01 (0.02)	0.01 (0.03)	1.000	-0.01 (0.05)	-0.01 (0.03)	-0.01 (0.04)	1.000
	Interstitial infiltrate classifier (i>1Prob)	-0.01 (0.04)	-0.02 (0.03)	-0.01 (0.04)	0.964	0.01 (0.04)	0.01 (0.04)	0.00 (0.04)	0.964	0.00 (0.05)	-0.01 (0.06)	-0.01 (0.05)	0.964
	Cytotoxic T cell-associated (QCAT)	-0.20 (0.41)	-0.70 (0.46)	-0.50 (0.43)	0.425	0.10 (0.33)	0.54 (0.57)	0.44 (0.45)	0.425	-0.10 (0.46)	-0.16 (0.41)	-0.06 (0.44)	0.871
	T cell burden (TCB)	-0.44 (0.71)	-0.79 (0.56)	-0.36 (0.64)	0.825	0.28 (0.57)	0.34 (0.6)	0.06 (0.58)	0.905	-0.16 (0.88)	-0.45 (0.43)	-0.29 (0.65)	0.825
macrophage-related	Alternatively activated macrophage (AMAT1)	0.08 (0.2)	-0.27 (0.45)	-0.35 (0.33)	0.128	0.09 (0.23)	0.11 (0.47)	0.02 (0.35)	0.812	0.18 (0.19)	-0.16 (0.32)	-0.33 (0.25)	0.128
	Constitutive macrophage (QCMAT)	0.00 (0.15)	-0.24 (0.27)	-0.24 (0.21)	0.340	0.01 (0.12)	0.09 (0.27)	0.08 (0.2)	0.685	0.01 (0.11)	-0.15 (0.24)	-0.16 (0.18)	0.348
injury-related	Injury-repair associated (IRRAT30)	0.18 (0.52)	0.05 (0.4)	-0.13 (0.46)	0.353	0.05 (0.66)	-0.11 (0.72)	-0.16 (0.69)	0.353	0.23 (0.6)	-0.06 (0.48)	-0.29 (0.54)	0.175
	Injury-repair induced, day 3 (IRITD3)	0.03 (0.19)	0.04 (0.12)	0.01 (0.15)	0.945	0.03 (0.18)	-0.18 (0.18)	-0.20 (0.18)	0.060	0.06 (0.19)	-0.13 (0.14)	-0.19 (0.16)	0.060
	Injury-repair induced, day 5 (IRITD5)	-0.01 (0.17)	0.04 (0.14)	0.05 (0.15)	0.932	0.07 (0.2)	-0.17 (0.21)	-0.24 (0.21)	0.093	0.06 (0.18)	-0.13 (0.17)	-0.18 (0.18)	0.093
	Fibrosis classifier (ci>1Prob)	-0.02 (0.29)	0.13 (0.32)	0.15 (0.3)	0.434	0.06 (0.25)	-0.08 (0.53)	-0.14 (0.39)	0.844	0.04 (0.25)	0.05 (0.37)	0.01 (0.31)	0.434
	Atrophy classifier (ct>1Prob)	0.01 (0.33)	0.12 (0.31)	0.12 (0.32)	0.587	0.09 (0.27)	-0.06 (0.44)	-0.15 (0.36)	0.587	0.10 (0.32)	0.07 (0.41)	-0.03 (0.36)	0.587
parenchyma-related	Kidney parenchymal (KT1)	-0.05 (0.22)	-0.15 (0.23)	-0.09 (0.22)	0.731	-0.01 (0.11)	0.10 (0.22)	0.11 (0.17)	0.705	-0.06 (0.25)	-0.04 (0.18)	0.02 (0.22)	0.705
	Kidney parenchymal - no solute carriers (KT2)	-0.09 (0.39)	-0.16 (0.37)	-0.07 (0.38)	0.577	-0.09 (0.2)	0.14 (0.38)	0.23 (0.29)	0.484	-0.18 (0.42)	-0.02 (0.31)	0.16 (0.37)	0.577
Grey shading denotes ANOVA interactive effect FDR < 0.05 FDR correction was carried out within each annotation grouping

1.7 PLOT THE RESULTS

# EXTRACT LEGEND FOR PLOTS ####
panel_legend <- plots_artANOVA %>%
    dplyr::filter(variable == "AMAT1") %>%
    pull(plot_violin) %>%
    ggpubr::get_legend() %>%
    ggpubr::as_ggplot() +
    theme(plot.margin = unit(c(0, 0, -1, 0), "cm"))


# MOLECULAR ABMR PANELS ####
panel_violin_abmr <- plots_artANOVA %>%
    dplyr::filter(category %in% c("ABMR")) %>%
    pull(plot_violin) %>%
    wrap_plots(nrow = 1, ncol = 5) &
    theme(
        axis.text = element_text(size = 10, colour = "black"),
        legend.position = "none",
        plot.background = element_rect(fill = "grey95", colour = " white")
    )

panel_violin_abmr <- panel_violin_abmr %>%
    ggarrange(
        labels = "B",
        font.label = list(size = 25, face = "bold"),
        legend = "none"
    ) %>%
    ggpubr::annotate_figure(
        top = text_grob("Effect of felzartamab on molecular ABMR activity scores",
            face = "bold.italic",
            size = 25,
            hjust = 1.27
        )
    )


# MOLECULAR TCMR PANELS ####
panel_violin_tcmr <- plots_artANOVA %>%
    dplyr::filter(category %in% c("TCMR")) %>%
    pull(plot_violin) %>%
    wrap_plots(nrow = 1, ncol = 5) &
    theme(
        axis.text = element_text(size = 10, colour = "black"),
        legend.position = "none",
    )

panel_violin_tcmr <- panel_violin_tcmr %>%
    ggarrange(
        labels = "C",
        font.label = list(size = 25, face = "bold"),
        legend = "none"
    ) %>%
    ggpubr::annotate_figure(
        top = text_grob("Effect of felzartamab on molecular TCMR activity scores",
            face = "bold.italic",
            size = 25,
            hjust = 1.27
        )
    )


# COMBINED VIOLIN PANELS ####
panels_violin <- ggarrange(
    panel_violin_abmr,
    panel_violin_tcmr,
    nrow = 2,
    heights = c(1, 1)
)

ggarrange(
    panel_legend,
    panels_violin,
    nrow = 2,
    heights = c(0.125, 1)
) 

2 Analysing the effect of treatment on gene expression

2.1 INTRODUCTION

This analysis assesses how felzartamab therapy affects gene expression from biopsies from kidney allografts taken from patients with antibody mediated rejection. We use the Linear Models for Microarray Data (limma) pipeline for applying empirical Bayes mediated t-tests on individual genes. Specific contrasts to assessed the interaction between treatment and time are specified to estimate the effect of treatemnt on each gene.

References:

• Ritchie, M.E., et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 43, e47 (2015).

2.2 DATA WRANGLING

First we need to load the necessary R packages and data.

# HOUSEKEEPING ####
# CRAN libraries
library(tidyverse) 
library(flextable) 
library(officer)
# Bioconductor libraries
library(Biobase) 
library(limma) 
library(genefilter) 
# Custom operators and functions
"%nin%" <- function(a, b) match(a, b, nomatch = 0) == 0
# Load data
load("C:/R/CD38-effect-of-treatment/natmed/data/cd38_3Sept24.RData")
load("C:/R/CD38-effect-of-treatment/natmed/data/affymap219.RData")

# Seed for reproducibility
set.seed(42)

# DEFINE THE SET ####
set00 <- set[, set$Patient %nin% c(15, 18)]

• We filter the gene expression data by to remove genes with minimal variance across the population (IQR filtering).

# IQR FILTER THE DATA ####
f1 <- function(x) (IQR(x) > 0.5)
ff <- filterfun(f1)
if (!exists("selected")) {
    selected <- genefilter(set00, ff)
}
set01 <- set00[selected, ]

• We then retain a single probeset for each gene based on whichever has the highest mean expression in the population.

# KEEP UNIQUE GENES (keep probe with highest mean expression) ####
mean_exprs_by_probe <- set01 %>%
    exprs() %>%
    as_tibble(rownames = "AffyID") %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb) %>% tibble(), ., by = "AffyID") %>%
    mutate(mean_exprs = set01 %>%
        exprs() %>% rowMeans(), .after = Symb)

genes <- mean_exprs_by_probe %>%
    group_by(Symb) %>%
    dplyr::slice_max(mean_exprs) %>%
    dplyr::filter(Symb != "") %>%
    distinct(Symb, .keep_all = TRUE) %>%
    pull(AffyID)

set <- set01[featureNames(set01) %in% genes, ]

2.3 REMOVE GENES THAT DIFFER AT BASELINE BETWEEN PLACEBO AND FELZARTAMAB ARMS

The felzartamab trial was a phase II randomize control trial (RCT). While RCTs are the ‘gold standard’ for establishing causality (i.e., the effect of treatment), they are not perfect, particularly when the number of study participants is low. Thus, we pre-screened gene expression data to remove any genes that differed between placebo and felzartamab arm as baseline. This was done to reduce the potential for interpreting the effect of treatment on genes with dissimilar baseline profiles.

# DEFINE FACTOR FOR CONTRASTS ####
Felzartamab_Followup <- set$Felzartamab_Followup %>% droplevels()


# BLOCK DESIGN week24 - baseline ####
design_block01 <- model.matrix(~ 0 + Felzartamab_Followup)
contrast_block_01 <- makeContrasts(
    "x =  (Felzartamab_FollowupBaseline_Felzartamab) -(Felzartamab_FollowupBaseline_Placebo)",
    levels = design_block01
)

# FIT BLOCK week24 - baseline LIMMA MODEL ####
fit_block_1 <- limma::lmFit(set, design_block01)
cfit_block_1 <- limma::contrasts.fit(fit_block_1, contrast_block_01)
ebayes_block_1 <- limma::eBayes(cfit_block_1)
tab_block_1 <- limma::topTable(ebayes_block_1, adjust = "fdr", sort.by = "p", number = "all")

# CALCULATE MEAN GENE EXPRESSION FOR EACH PROBE BETWEEN GROUPINGS ####
means_baseline <- fit_block_1 %>%
    avearrays() %>%
    data.frame() %>%
    rownames_to_column("AffyID") %>%
    tibble() %>%
    mutate_if(is.numeric, ~ 2^. %>% round(0)) %>%
    rename_at(vars(contains("Felz")), ~ str_remove(., "Felzartamab_Followup")) %>%
    dplyr::select(AffyID, contains("Baseline"))

The top 20 genes that differ at baseline (sorted by p-value) between placebo and felzartamab arms are summarized below.

(Note that all FDR > 0.05. This is a consequence of differential expression among with n = 10 by treatment group across a large selection of genes (5,267 in this case). Nonetheless, genes with p < 0.05 will be excluded from subsequent analyses).

# FORMAT TOPTABLES ####
baseline_deg <- tab_block_1 %>%
    as_tibble(rownames = "AffyID") %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID") %>%
    arrange(P.Value) %>%
    mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
    tibble() %>%
    left_join(means_baseline, by = "AffyID") %>%
    dplyr::select(
        AffyID, Symb,
        Baseline_Placebo, Baseline_Felzartamab,
        t, logFC, P.Value, adj.P.Val,
    ) %>%
    dplyr::rename(p = P.Value, FDR = adj.P.Val)

baseline_deg %>%
    dplyr::slice_min(p, n = 20) %>%
    mutate_at(vars(t, logFC, FDR), ~ formatC(.x, format = "f", digits = 2)) %>%
    mutate_at(vars(p), ~ formatC(.x, format = "f", digits = 5)) %>%
    flextable::flextable()

AffyID	Symb	Baseline_Placebo	Baseline_Felzartamab	t	logFC	p	FDR
11727331_at	SPAG4	90	50	-4.23	-0.85	0.00008	0.33
11763725_a_at	HNRNPU	366	665	4.10	0.86	0.00013	0.33
11718286_a_at	UBL3	74	142	3.66	0.94	0.00053	0.45
11732590_x_at	TRIM5	42	63	3.65	0.59	0.00055	0.45
11725671_a_at	IPO7	42	66	3.54	0.64	0.00079	0.45
11741431_a_at	KCNJ16	203	311	3.49	0.61	0.00091	0.45
11726915_at	SLCO4C1	87	176	3.48	1.02	0.00096	0.45
11723109_a_at	GPBP1	153	263	3.46	0.78	0.00099	0.45
11739744_x_at	MXI1	123	191	3.45	0.64	0.00102	0.45
11720890_a_at	AGL	39	63	3.40	0.70	0.00121	0.45
11746431_a_at	FBXO11	28	45	3.38	0.68	0.00127	0.45
11751841_a_at	TNFRSF17	56	20	-3.36	-1.51	0.00134	0.45
11744367_a_at	ZMYND11	30	49	3.33	0.68	0.00150	0.45
11748857_a_at	AK3	371	578	3.32	0.64	0.00155	0.45
11727307_x_at	ZNF320	46	66	3.31	0.51	0.00160	0.45
11747648_a_at	TCAIM	96	152	3.29	0.66	0.00168	0.45
11728966_s_at	ZNF28	29	67	3.27	1.21	0.00181	0.45
11723387_a_at	PTAR1	149	261	3.25	0.81	0.00187	0.45
11739771_a_at	OXR1	31	51	3.24	0.73	0.00194	0.45
11730038_at	GNG13	44	28	-3.24	-0.63	0.00196	0.45

2.4 ESTIMATING TREATMENT EFFECT ON GENE EXPRESSION

We now carry out differential expression analysis to estimate the treatment effect on individual genes.

# DEFINE GENES SIMILAR AT BASELINE ####
genes_baseline <- baseline_deg %>%
    dplyr::filter(p > 0.05) %>%
    pull(AffyID)

# WRANGLE THE MEAN EXPRESSION DATA ####
mean_exprs_by_probe <- set01 %>%
    exprs() %>%
    as_tibble(rownames = "AffyID") %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb) %>% tibble(), ., by = "AffyID") %>%
    mutate(
        mean_exprs = set01 %>%
            exprs() %>%
            rowMeans(),
        .after = Symb
    )

genes <- mean_exprs_by_probe %>%
    group_by(Symb) %>%
    dplyr::slice_max(mean_exprs) %>%
    dplyr::filter(Symb != "", AffyID %in% genes_baseline) %>%
    distinct(Symb, .keep_all = TRUE) %>%
    pull(AffyID)

set02 <- set01[featureNames(set01) %in% genes, ]

We need to define the model structure for the contrasts we intend to make. Here we will define 3 major contrast groups; 1) effect of time in placebo group, 2) effect of time in felzartamab group, and 3) the interactive effect of time and treatment. We do this for each time window (baseline to week 24, week 24 to week 52, and baseline to week 52)

# DEFINE FACTOR FOR CONTRASTS ####
Felzartamab_Followup <- set02$Felzartamab_Followup %>% droplevels()
cortex <- set02$Cortexprob


# DESIGN ####
design <- model.matrix(~ 0 + Felzartamab_Followup + cortex)


# CONTRAST DESIGN week24 - baseline ####
contrast_interaction_01 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek24_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2 -(Felzartamab_FollowupWeek24_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
    levels = design
)
contrast_felzartamab_01 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek24_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2",
    levels = design
)
contrast_placebo_01 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek24_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
    levels = design
)


# CONTRAST DESIGN week52 - week24 ####
contrast_interaction_02 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupWeek24_Felzartamab)/2 - (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupWeek24_Placebo)/2",
    levels = design
)
contrast_felzartamab_02 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupWeek24_Felzartamab)/2",
    levels = design
)
contrast_placebo_02 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupWeek24_Placebo)/2",
    levels = design
)


# CONTRAST DESIGN week52 - baseline####
contrast_interaction_03 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2 - (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
    levels = design
)
contrast_felzartamab_03 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2",
    levels = design
)
contrast_placebo_03 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
    levels = design
)

We now fit the models we’ve specified

# FIT BLOCK week24 - baseline LIMMA MODEL ####
fit_interaction_1 <- limma::lmFit(set02, design)
cfit_interaction_1 <- limma::contrasts.fit(fit_interaction_1, contrast_interaction_01)
ebayes_interaction_1 <- limma::eBayes(cfit_interaction_1)
tab_interaction_1 <- limma::topTable(ebayes_interaction_1, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_interaction_1$stdev.unscaled * sqrt(ebayes_interaction_1$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")

fit_felzartamab_1 <- limma::lmFit(set02, design)
cfit_felzartamab_1 <- limma::contrasts.fit(fit_felzartamab_1, contrast_felzartamab_01)
ebayes_felzartamab_1 <- limma::eBayes(cfit_felzartamab_1)
tab_felzartamab_1 <- limma::topTable(ebayes_felzartamab_1, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_felzartamab_1$stdev.unscaled * sqrt(ebayes_felzartamab_1$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")

fit_placebo_1 <- limma::lmFit(set02, design)
cfit_placebo_1 <- limma::contrasts.fit(fit_placebo_1, contrast_placebo_01)
ebayes_placebo_1 <- limma::eBayes(cfit_placebo_1)
tab_placebo_1 <- limma::topTable(ebayes_placebo_1, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_placebo_1$stdev.unscaled * sqrt(ebayes_placebo_1$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")


# FIT BLOCK week52 - week24 LIMMA MODEL ####
fit_interaction_2 <- limma::lmFit(set02, design)
cfit_interaction_2 <- limma::contrasts.fit(fit_interaction_2, contrast_interaction_02)
ebayes_interaction_2 <- limma::eBayes(cfit_interaction_2)
tab_interaction_2 <- limma::topTable(ebayes_interaction_2, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_interaction_2$stdev.unscaled * sqrt(ebayes_interaction_2$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")

fit_felzartamab_2 <- limma::lmFit(set02, design)
cfit_felzartamab_2 <- limma::contrasts.fit(fit_felzartamab_2, contrast_felzartamab_02)
ebayes_felzartamab_2 <- limma::eBayes(cfit_felzartamab_2)
tab_felzartamab_2 <- limma::topTable(ebayes_felzartamab_2, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_felzartamab_2$stdev.unscaled * sqrt(ebayes_felzartamab_2$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")

fit_placebo_2 <- limma::lmFit(set02, design)
cfit_placebo_2 <- limma::contrasts.fit(fit_placebo_2, contrast_placebo_02)
ebayes_placebo_2 <- limma::eBayes(cfit_placebo_2)
tab_placebo_2 <- limma::topTable(ebayes_placebo_2, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_placebo_2$stdev.unscaled * sqrt(ebayes_placebo_2$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")


# FIT BLOCK week52 - baseline LIMMA MODEL ####
fit_interaction_3 <- limma::lmFit(set02, design)
cfit_interaction_3 <- limma::contrasts.fit(fit_interaction_3, contrast_interaction_03)
ebayes_interaction_3 <- limma::eBayes(cfit_interaction_3)
tab_interaction_3 <- limma::topTable(ebayes_interaction_3, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_interaction_3$stdev.unscaled * sqrt(ebayes_interaction_3$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")

fit_felzartamab_3 <- limma::lmFit(set02, design)
cfit_felzartamab_3 <- limma::contrasts.fit(fit_felzartamab_3, contrast_felzartamab_03)
ebayes_felzartamab_3 <- limma::eBayes(cfit_felzartamab_3)
tab_felzartamab_3 <- limma::topTable(ebayes_felzartamab_3, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_felzartamab_3$stdev.unscaled * sqrt(ebayes_felzartamab_3$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")

fit_placebo_3 <- limma::lmFit(set02, design)
cfit_placebo_3 <- limma::contrasts.fit(fit_placebo_3, contrast_placebo_03)
ebayes_placebo_3 <- limma::eBayes(cfit_placebo_3)
tab_placebo_3 <- limma::topTable(ebayes_placebo_3, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_placebo_3$stdev.unscaled * sqrt(ebayes_placebo_3$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")

2.5 ORAGNIZE THE RESULTS

The results are now organized for tabulation

# MERGE BLOCKS ####
tab_block_1 <- tab_interaction_1 %>%
    left_join(
        tab_felzartamab_1 %>%
            dplyr::select(AffyID, logFC) %>%
            rename(flogFC = logFC),
        by = "AffyID"
    ) %>%
    left_join(
        tab_placebo_1 %>%
            dplyr::select(AffyID, logFC) %>%
            rename(plogFC = logFC),
        by = "AffyID"
    ) %>%
    relocate(plogFC, flogFC, .before = logFC) %>%
    arrange(P.Value)

tab_block_2 <- tab_interaction_2 %>%
    left_join(
        tab_felzartamab_2 %>%
            dplyr::select(AffyID, logFC) %>%
            rename(flogFC = logFC),
        by = "AffyID"
    ) %>%
    left_join(
        tab_placebo_2 %>%
            dplyr::select(AffyID, logFC) %>%
            rename(plogFC = logFC),
        by = "AffyID"
    ) %>%
    relocate(plogFC, flogFC, .before = logFC) %>%
    arrange(P.Value)

tab_block_3 <- tab_interaction_3 %>%
    left_join(
        tab_felzartamab_3 %>%
            dplyr::select(AffyID, logFC) %>%
            rename(flogFC = logFC),
        by = "AffyID"
    ) %>%
    left_join(
        tab_placebo_3 %>%
            dplyr::select(AffyID, logFC) %>%
            rename(plogFC = logFC),
        by = "AffyID"
    ) %>%
    relocate(plogFC, flogFC, .before = logFC) %>%
    arrange(P.Value)


# CALCULATE MEAN GENE EXPRESSION FOR EACH PROBE BETWEEN GROUPINGS ####
means <- fit_interaction_1 %>%
    avearrays() %>%
    as_tibble(rownames = "AffyID") %>%
    mutate_if(is.numeric, ~ 2^. %>% round(0)) %>%
    rename_at(vars(contains("Felz")), ~ str_remove(., "Felzartamab_Followup")) %>%
    dplyr::select(-contains("Week12"), -any_of(c("cortex")))


# FORMAT TOPTABLES ####
table_interaction_1 <- tab_block_1 %>%
    arrange(P.Value) %>%
    mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
    tibble() %>%
    left_join(means, by = "AffyID") %>%
    dplyr::select(
        AffyID, Symb, Gene, PBT,
        t, plogFC, flogFC, logFC, se, P.Value, adj.P.Val,
        all_of(colnames(means))
    ) %>%
    dplyr::rename(
        p = P.Value,
        FDR = adj.P.Val
    )

table_interaction_2 <- tab_block_2 %>%
    arrange(P.Value) %>%
    mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
    tibble() %>%
    left_join(means, by = "AffyID") %>%
    dplyr::select(
        AffyID, Symb, Gene, PBT,
        t, plogFC, flogFC, logFC, se, P.Value, adj.P.Val,
        all_of(colnames(means))
    ) %>%
    dplyr::rename(
        p = P.Value,
        FDR = adj.P.Val
    )

table_interaction_3 <- tab_block_3 %>%
    arrange(P.Value) %>%
    mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
    tibble() %>%
    left_join(means, by = "AffyID") %>%
    dplyr::select(
        AffyID, Symb, Gene, PBT,
        t, plogFC, flogFC, logFC, se, P.Value, adj.P.Val,
        all_of(colnames(means))
    ) %>%
    dplyr::rename(
        p = P.Value,
        FDR = adj.P.Val
    )

deg <- tibble(
    design = c(
        "Baseline_vs_Week24",
        "Week24_vs_Week52",
        "Baseline_vs_Week52"
    ),
    table = list(
        table_interaction_1,
        table_interaction_2,
        table_interaction_3
    )
)

2.6 SUMMARIZE THE RESULTS

# FORMAT TABLES TO MAKE FLEXTABLES ####
table_deg <- deg %>%
    expand_grid(direction = c("all", "increased", "decreased")) %>%
    relocate(direction, .after = design) %>%
    mutate(
        gene_tables = pmap(
            list(direction, table),
            function(direction, table) {
                if (direction == "increased") {
                    table <- table %>%
                        dplyr::filter(logFC > 0) %>%
                        arrange(p)
                } else if (direction == "decreased") {
                    table <- table %>%
                        dplyr::filter(logFC < 0) %>%
                        arrange(p)
                }
                table %>%
                    dplyr::select(-t, -se, -FDR, -contains("AffyID"), -contains("MMDx")) %>%
                    dplyr::slice(1:20) %>%
                    mutate(
                        Gene = Gene %>% str_remove("///.*"),
                        plogFC = plogFC %>% round(2),
                        flogFC = flogFC %>% round(2),
                        logFC = logFC %>% round(2),
                        p = case_when(
                            p < 0.0001 ~ p %>% formatC(digits = 0, format = "e"),
                            TRUE ~ p %>% formatC(digits = 4, format = "f")
                        )
                    )
            }
        )
    )


# GLOBAL PARAMETERS FOR FLEXTABLES ####
header1 <- c(
    "Gene\nsymbol", "Gene", "PBT",
    rep("Differential expression", 4),
    rep("Mean expression by group", 6)
)
header2 <- c(
    "Gene\nsymbol", "Gene", "PBT",
    "\u394\nplacebo\nlogFC", "\u394\nfelzartamab\nlogFC",
    "\u394\u394\nlogFC", "\u394\u394\nP",
    rep("Placebo", 3), rep("Felzartamab", 3)
)
header3 <- c(
    "Gene\nsymbol", "Gene", "PBT",
    "\u394\nplacebo\nlogFC", "\u394\nfelzartamab\nlogFC",
    "\u394\u394\nlogFC", "\u394\u394\nP",
    "Baseline\n(N=10)", "Week24\n(N=10)", "Week52\n(N=10)", "Baseline\n(N=10)", "Week24\n(N=10)", "Week52\n(N=10)"
)

# MAKE FORMATTED FLEXTABLES ####
flextables <- table_deg %>%
    mutate(
        flextable = pmap(
            list(design, gene_tables),
            function(design, gene_tables) {
                if (design == "Baseline_vs_Week24") {
                    title <- paste("Table i. Top 20 differentially expressed genes between baseline and week24 in biopsies from placebo and Felzartamab treated patients (by P-value)", sep = "")
                } else if (design == "Week24_vs_Week52") {
                    title <- paste("Table i. Top 20 differentially expressed genes between week24 and week52 in biopsies from placebo and Felzartamab treated patients (by P-value)", sep = "")
                } else if (design == "Baseline_vs_Week52") {
                    title <- paste("Table i. Top 20 differentially expressed genes between baseline and week52 in biopsies from placebo and Felzartamab treated patients (by P-value)", sep = "")
                }
                gene_tables %>%
                    arrange(logFC) %>%
                    flextable::flextable() %>%
                    flextable::delete_part("header") %>%
                    flextable::add_header_row(top = TRUE, values = header3) %>%
                    flextable::add_header_row(top = TRUE, values = header2) %>%
                    flextable::add_header_row(top = TRUE, values = header1) %>%
                    flextable::add_header_row(top = TRUE, values = rep(title, ncol_keys(.))) %>%
                    flextable::merge_v(part = "header") %>%
                    flextable::merge_h(part = "header") %>%
                    flextable::border_remove() %>%
                    flextable::border(part = "header", border = officer::fp_border()) %>%
                    flextable::border(part = "body", border = officer::fp_border()) %>%
                    flextable::align(align = "center") %>%
                    flextable::align(align = "center", part = "header") %>%
                    flextable::font(fontname = "Arial", part = "all") %>%
                    flextable::fontsize(size = 7, part = "all") %>%
                    flextable::fontsize(size = 7, part = "footer") %>%
                    flextable::fontsize(i = 1, size = 12, part = "header") %>%
                    flextable::bold(part = "header") %>%
                    flextable::bg(bg = "white", part = "all") %>%
                    flextable::padding(padding = 0, part = "all") %>%
                    flextable::autofit()
            }
        )
    )

Genes suppressed by felzartamab in week 24 compared to baseline
(Note that all ΔΔFDR > 0.05. This is a consequence of differential expression among with n = 10 by treatment group across a large selection of genes (5,267 in this case). Thus, treatment effects on individual genes should not be interpreted. Instead, geneset enrichment analysis (GSEA) will be used to assess bulk trends).

flextables %>%
    dplyr::filter(design == "Baseline_vs_Week24", direction == "decreased") %>%
    pull(flextable) %>%
    pluck(1)

Table i. Top 20 differentially expressed genes between baseline and week24 in biopsies from placebo and Felzartamab treated patients (by P-value)
Gene symbol	Gene	PBT	Differential expression				Mean expression by group
			Δ placebo logFC	Δ felzartamab logFC	ΔΔ logFC	ΔΔ P	Placebo			Felzartamab
							Baseline (N=10)	Week24 (N=10)	Week52 (N=10)	Baseline (N=10)	Week24 (N=10)	Week52 (N=10)
CXCL11	chemokine (C-X-C motif) ligand 11	ABMR-RAT,GRIT3,RAT,Rej-RAT	-0.16	-1.01	-0.85	0.0176	494	395	372	617	153	317
CXCL9	chemokine (C-X-C motif) ligand 9	ABMR-RAT,GRIT1,GRIT3,RAT,Rej-RAT	-0.14	-0.94	-0.79	0.0290	1,344	1,102	1,384	1,684	459	1,036
CXCL10	chemokine (C-X-C motif) ligand 10	ABMR-RAT,GRIT1,GRIT3,RAT,Rej-RAT	-0.21	-0.98	-0.77	0.0228	844	632	720	1,063	274	615
IDO1	indoleamine 2,3-dioxygenase 1	ABMR-RAT,GRIT3,RAT,Rej-RAT	-0.09	-0.77	-0.68	0.0103	338	298	296	409	140	251
LYZ	lysozyme	IRITD5,QCMAT	-0.01	-0.67	-0.66	0.0386	1,676	1,658	1,884	2,244	884	1,184
FCGR3A	Fc fragment of IgG, low affinity IIIa, receptor (CD16a)	ABMR-RAT,GRIT2,IRRAT950,RAT,Rej-RAT	0.06	-0.58	-0.64	0.0167	550	602	613	571	256	352
GZMB	granzyme B	ABMR-RAT,QCAT,RAT,Rej-RAT	-0.01	-0.62	-0.61	0.0019	149	147	130	143	61	114
GBP1	guanylate binding protein 1, interferon-inducible	ABMR-RAT,GRIT3,RAT,Rej-RAT	-0.03	-0.61	-0.59	0.0096	678	655	643	909	389	556
GNLY	granulysin	ABMR-RAT,DSAST,QCAT,RAT,Rej-RAT	-0.18	-0.75	-0.57	0.0285	192	149	134	165	58	100
SH2D1B	SH2 domain containing 1B	ABMR-RAT,DSAST,NKB,RAT	0.07	-0.44	-0.51	0.0010	38	42	34	40	22	33
FCGR1A	Fc fragment of IgG, high affinity Ia, receptor (CD64)	GRIT2,IRRAT950,RAT,TCMR-RAT	0.09	-0.41	-0.50	0.0304	131	148	161	152	86	86
KLRD1	killer cell lectin-like receptor subfamily D, member 1	ABMR-RAT,RAT,Rej-RAT	-0.06	-0.49	-0.43	0.0295	110	101	96	126	63	102
EPB41L3	erythrocyte membrane protein band 4.1-like 3		0.22	-0.20	-0.42	0.0388	209	286	258	339	258	283
LOC339059	uncharacterized LOC339059		0.26	-0.15	-0.41	0.0203	51	73	55	66	53	59
APOL1	apolipoprotein L1	ABMR-RAT,GRIT3,RAT,Rej-RAT	0.07	-0.32	-0.39	0.0205	626	687	744	754	483	636
LYPD5	LY6/PLAUR domain containing 5	ABMR-RAT,RAT,Rej-RAT	0.04	-0.33	-0.37	0.0255	27	28	27	30	19	26
TRAV12-2	T cell receptor alpha variable 12-2		0.25	-0.10	-0.35	0.0173	39	55	39	49	43	43
IL18RAP	interleukin 18 receptor accessory protein		0.02	-0.32	-0.34	0.0233	38	39	41	42	27	33
PRF1	perforin 1 (pore forming protein)	ABMR-RAT,QCAT,RAT,Rej-RAT	-0.02	-0.35	-0.33	0.0382	152	148	137	145	89	123
PMM2	phosphomannomutase 2		0.20	-0.07	-0.27	0.0380	38	50	41	46	42	49

Genes suppressed by felzartamab in week 52 compared to baseline
(Note that all ΔΔFDR > 0.05. This is a consequence of differential expression among with n = 10 by treatment group across a large selection of genes (5,267 in this case). Thus, treatment effects on individual genes should not be interpreted. Instead, geneset enrichment analysis (GSEA) will be used to assess bulk trends).

flextables %>%
    dplyr::filter(design == "Baseline_vs_Week52", direction == "decreased") %>%
    pull(flextable) %>%
    pluck(1)

Table i. Top 20 differentially expressed genes between baseline and week52 in biopsies from placebo and Felzartamab treated patients (by P-value)
Gene symbol	Gene	PBT	Differential expression				Mean expression by group
			Δ placebo logFC	Δ felzartamab logFC	ΔΔ logFC	ΔΔ P	Placebo			Felzartamab
							Baseline (N=10)	Week24 (N=10)	Week52 (N=10)	Baseline (N=10)	Week24 (N=10)	Week52 (N=10)
HNRNPK	heterogeneous nuclear ribonucleoprotein K		0.28	-0.38	-0.66	0.0055	818	1,044	1,203	1,096	1,059	649
COL1A1	collagen, type I, alpha 1	COLA1356,FICOL,IRITD5,LivGST_UP	0.38	-0.26	-0.64	0.0075	310	333	526	329	320	230
NFIX	nuclear factor I/X (CCAAT-binding transcription factor)	CT1	0.42	-0.15	-0.57	0.0012	162	179	289	221	167	180
RHOA	ras homolog family member A	IRRAT950	0.28	-0.29	-0.57	0.0049	691	803	1,017	923	809	620
COL1A2	collagen, type I, alpha 2	FICOL,IRITD5,LivGST_UP	0.26	-0.31	-0.57	0.0118	1,238	1,152	1,769	1,534	1,194	997
ITGB3	integrin beta 3	IRRAT30,IRRAT950	0.34	-0.21	-0.55	0.0100	105	137	169	110	103	82
LOC100129518	uncharacterized LOC100129518	GRIT3,IRRAT950	0.25	-0.22	-0.47	0.0050	429	546	605	527	461	390
ACTB	actin, beta	IRRAT950	0.15	-0.31	-0.47	0.0077	6,451	7,352	7,993	7,274	5,662	4,720
SPARC	secreted protein, acidic, cysteine-rich (osteonectin)	IRITD3	0.16	-0.30	-0.46	0.0100	2,411	2,532	3,011	2,343	2,094	1,550
DNAJA4	DnaJ (Hsp40) homolog, subfamily A, member 4	IRRAT950	0.27	-0.14	-0.42	0.0095	49	68	72	68	61	56
SLFN5	schlafen family member 5		0.19	-0.20	-0.40	0.0050	39	44	51	49	38	37
CYR61	cysteine-rich, angiogenic inducer, 61	IRITD3,IRRAT950	0.14	-0.26	-0.40	0.0074	271	268	330	338	234	236
SDCBP	syndecan binding protein		0.24	-0.14	-0.38	0.0120	1,135	1,404	1,582	1,478	1,358	1,224
SBSPON	somatomedin B and thrombospondin type 1 domain containing		0.31	-0.08	-0.38	0.0123	217	250	332	188	249	169
PAFAH1B2	platelet-activating factor acetylhydrolase 1b, catalytic subunit 2 (30kDa)		0.23	-0.13	-0.36	0.0099	309	372	423	386	396	323
TMEM2	transmembrane protein 2		0.10	-0.26	-0.36	0.0118	270	307	311	318	259	222
CALM2	calmodulin 2 (phosphorylase kinase, delta)	CT1,IRITD3	0.12	-0.25	-0.36	0.0138	1,304	1,427	1,538	1,623	1,342	1,155
WNK1	WNK lysine deficient protein kinase 1		0.26	-0.08	-0.34	0.0104	703	884	1,012	894	996	799
CEBPB	CCAAT/enhancer binding protein (C/EBP), beta	IRITD3	0.11	-0.22	-0.33	0.0123	419	452	488	417	333	307
STAT6	signal transducer and activator of transcription 6, interleukin-4 induced		0.09	-0.20	-0.30	0.0144	671	710	765	746	616	563